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Abstract 

We discuss analytical approximation schemes for the dynamics of diluted spin models. 
The original dynamics of the complete set of degrees of freedom is replaced by a hierarchy of 
equations including an increasing number of global observables, which can be closed approxi- 
mately at different levels of the hierarchy. We illustrate this method on the simple example of 
the Ising ferromagnet on a Bethe lattice, investigating the first three possible closures, which 
are all exact in the long time limit, and which yield more and more accurate predictions for 
the finite-time behavior. We also investigate the critical region around the phase transition, 
and the behavior of two-time correlation functions. We finally underline the close relation- 
ship between this approach and the dynamical replica theory under the assumption of replica 
symmetry. 

LPT-ENS 04/04 

1 Introduction 

The last few years have seen a considerable increase in the research activity concerning diluted and 
disordered spin models. This recent interest is mainly driven by two very different, but technically 
closely related motivations. 

The first one aims at understanding common fundamental features of glassy systems, or more 
generally of out-of-equilibrium problems. Besides various other approaches ^ the study of 
disordered spin models has led to important insights and given rise to the use of general concepts 
like, e.g., replica symmetry breaking on the static side or effective temperatures on the dynamic 
one 0J|S]. Most of these studies were performed on mean-field models, and the validity of these 
concepts for finite-dimensional systems is still under discussion. For spin models, mean-field has 
long been another name for fully-connected, i.e. models where each degree of freedom interacts 
with each other. More recently, a lot of effort was invested in the study of diluted models El 
|H| El EH 1111 El • These are still mean-field like in the sense that they do not have an underlying 
geometric structure, and are thus analytically easier to treat than finite-dimensional problems. At 
the same time, they have finite connectivity: each degree of freedom interacts only with a finite 
number of neighbors, and the concept of a local environment is well-defined. This allows, e.g., to 
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introduce microscopically motivated interactions. Diluted models can therefore be considered as 
an intermediate step between the fully-connected and the finite-dimensional case. 

A second major motivation for the study of diluted spin models arises from their close connec- 
tion with a very interesting class of combinatorial optimization problems |13|. including famous 
problems like satisfiability of Boolean formulas and graph coloring for instance. These problems 
show phase transitions in the statistical properties of their solutions as well as in the dynamical 
behavior of algorithms and have thus been extensively studied using tools and concepts 

originally devised for statistical mechanics purpose \H\ El El 1201 IZH 122] ■ 

Despite the large interest in diluted models, many basic questions are still open. Whereas the 
main technical obstacles in analyzing the static behavior are solved by now, and the equilibrium 
properties of these models are quite well-understood, the knowledge on the dynamical side [231 
1241 1251 123 is still relatively poor. Here we reinvestigate some approximative methods recently 
introduced for analyzing the behavior of stochastic local search optimization algorithms [23 12H1 • 
We apply and generalize these ideas to a purely physical model, in order to better understand 
the assumptions behind the approximations made, and to assess their quality. We concentrate 
on a very simple system, namely a ferromagnetic Ising model defined on a diluted network of 
fixed connectivity, or Bethe lattice. This example allows to easily present, test and understand 
approximate methods which should be useful to treat more interesting, glassy problems. In fact the 
non-equilibrium flavor of this work comes from the study of transient relaxation from an arbitrary 
initial configuration to thermal equilibrium, and not from the absence of thermal equilibrium as 
in the case of glassy systems, or from the lack of detailed balance as in previous investigations of 
algorithms. 

The paper is organized as follows. After this general introduction, the model and its equilibrium 
behavior are presented in Sec. [21 In Sec.|31the dynamical rules of the model are presented, and a 
hierarchical approximation scheme is developed. The first three approximations are compared to 
numerical simulations in Sec. ^ and they are further exploited analytically to describe the critical 
behavior in Sec. Sec. is dedicated to an extension of the previous approximations to the 
analysis of two-time quantities. In Sec. [3 we unveil the relation of our approach to the dynamical 
replica theory proposed by Coolen and Sherrington [23 CHI • Finally, the results are summarized 
in the last section, and possible future directions of research are presented. 

2 The model and its equilibrium behavior 

We consider a ferromagnetic Ising model on a Bethe lattice, given by its Hamiltonian 



depending on the microscopic configuration a = (<o~\, <jn) of the N Ising variables <7j = ±1, i = 
1, N. The couplings Jtj take the value +1 whenever two spin <Ji and ctj are connected, and zero 
otherwise. Compared to the usual form we have shifted the Hamiltonian by its ground state value, 
and divided it by two. In this way, the Hamiltonian simply counts the number of unsatisfied edges, 
i.e. the number of edges carrying anti-parallel spins on their extremities. It can thus be rewritten 
as 



The reason why we choose this slightly modified form will become clear below, it allows a simpler 
presentation of the dynamical equations. 

As already said above, we consider this Ising model on a Bethe lattice. In accordance with the 
recent use in statistical mechanics 0, we define these as random regular graphs, i.e. all vertices 
have the same degree (number of neighboring sites) E21 ■ These graphs are locally tree- like. 
In contrast to Cayley trees they have no boundary, all vertices have exactly K neighbors and are 
thus equivalent. The graph therefore contains loops for any K > 1. These are, however, of length 
OQnN), i.e. they become long in the thermodynamic limit N — > oo (with K kept constant). Note 
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that the randomness of these graphs appears only through these loops. On finite length scales they 
appear homogeneous due to their constant vertex degree. 

Before investigating the non-equilibrium dynamics of this model, we shortly review its static 
behavior which can be easily solved using the Bethe-Peierls iterative approach |331 17] . For doing 
so, we assume for a moment that our model is defined on a tree, the long loops will be taken into 
account later as self-consistency conditions. 

Consider a given bond i.e. Jy = 1. Let us denote Z^icTi) the partition function of the 

subtree rooted in i, with deleted, and with a fixed value of spin <7j. These partition functions 
can be easily computed iteratively, 



Zi\j(ai)= Y{ ( E Z k\i( a k) exp{-/3 5 ai _ CTfe } J 
kjtj\J ik =i \o-k=±i / 



(3) 



where 8 denotes the inverse temperature. Defining the cavity field 

1 /%(+!) 
23 \Z lU (-l) 

we obtain the iterative equations 



2 ^ = E ln ( e .(- 1+ ^) +e -^„ ) ■ (b) 



At this point, we take into account the fact that the model is not defined on a tree, but that 
all vertices have the same vertex degree and are thus equivalent. We are therefore looking for a 
self-consistent homogeneous solution of Eq. (J5J with h = h^ij for all bonds 

, K — 1 ( e? h + e-W+fc) \ 

■ ln TrTZTTT, , —Ri ■ ( 6 ) 



2/3 \ e P(-i+h) +e -Ph 

This equation has the obvious paramagnetic solution h — 0. In fact, at high temperature this 
solution is the only one. There appears, however, a ferromagnetic phase transition at the critical 
inverse temperature /3 C = \n(K/(K — 2)). For lower temperatures, the trivial solution is ther- 
modynamically unstable, and two equivalent ferromagnetic solutions ±h describe the equilibrium 
behavior of the model. 

Here and in the following, without any loss of generality, we concentrate only on non-negative 
h. Once its value is known, we can immediately compute the Bethe free-energy density (total 
free-energy divided by N), 

13 f = (K-1) ln ^2 cosh (PjfZjh) ) ~ f" ln ( 2 C0S H 2 P h ) + . (7) 

Other interesting observables are the energy density 

d(3 \ K -1 ) 2 cosh(2(3h) + e-P w 

which equals 1^(1 + e* 9 ) -1 in the paramagnetic phase, and the magnetization 

m = tanh (b ^ ^ hj , (9) 

which becomes positive for positive cavity fields h. Note that the magnetization m depends on the 
cavity field h via the "true" effective field Kh/(K — 1) which incorporates the contributions of all 
K neighbors of a spin, cf. Eqs. ijoKill . 
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To explain some steps of the dynamical approach, we also need the probability p a {u) that a 
randomly selected vertex has spin value a and belongs to exactly u unsatisfied edges, i.e. u out of 
his K neighbors have spin —a. This quantity is given at equilibrium by 



p_(«) = ±r ( K \ e -^-Wh{K-u) ^ (1Q) 



7V„ 



with A4 being a normalization constant enforcing ^ M a p a (u) = 1. This can be rewritten in terms 
of the energy density and the magnetization as 

p»= 1 ^(f)( (1+ 2 ; m)if ) ( i - (1+ 2 ^) ■ (n) 

Similar expressions can be easily derived for the probability Pa 1 a 2 u 2) that a randomly selected 
edge has a first (resp. second) end vertex of spin o\ (resp. 02), and this vertex belongs to u\ (resp. 
U2) unsatisfied edges, 

A/ e \«i - 1/ W - 1/ 
p_+(ui,ua) = 1 ^ " 1 - !> ) e - M u 1+ » 2 -i)-2^(Jf-«+™-D 

jVe \tti - 1/ V"2 - IJ 

P—faUi) = J_^- 1 U if - 1 V^(« 1 +«2)-2 W 2K-„ 1 - U2 -2) (12) 
Me V "I / V U 2 / 

where M e is again determined by normalization. 



3 Dynamical approximation schemes 
3.1 Definitions 

We will study the following local stochastic dynamics of the model: in each algorithmic step 
T — > T + 1, a site i is chosen at random, i G {1, N}. This site is characterized by its spin Ui 
and the number Ui(a) = J"^- JijS au — aj of its unsatisfied incident edges. The spin is flipped to —oi 
with probability W(«i(<?), 0). We denote the new configuration, with spin <n flipped, by F%a. 

Obviously all edges incident to site i which were unsatisfied become satisfied, and vice versa. 
Hence the variation of the energy if the spin is flipped equals AE = K — 2ui . In order to reach 
thermal equilibrium at inverse temperature (3 in the long-time limit, we impose the detailed balance 
condition under the form 



W(u, (3) = W{K - u, 0) exp(-f3(K - 2u)) . (13) 
The two best-known possibilities falling in this category are the Metropolis rate, 

W(u,j3) =Tmn(l,e-^ K - 2u ^ , (14) 

and the Glauber rate, 

W(u,(3) = Ul-t & nh( (3 (^-u))) . (15) 



We will keep, however, a general form for W in our analytical treatment, only assuming Eq. (|13|) 
to be valid. In the thermodynamic limit, this discrete process acquires a continuous form by 
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defining the time as t = T/N, and stepwise differences of extensive observabies transiate into time 
derivatives of the corresponding observable densities (extensive observabies divided by N) . 

The dynamics of the system, or more precisely of the probabilities V(a,t) that a microscopic 
configuration a is found at time t, can be completely described by the master equation 

d 1 N 

-V(a, t) = ^J2 [-^M ^ 0) TO, *) + W( Ui (Fia), (3) V{F0, t)} . (16) 

i—l 

It is, however, far too complicated to solve these 2 N coupled equations directly. We are therefore 
going to present several approximative characterizations of the dynamics of the model. They all rely 
on the same idea: instead of following the evolution of the full distribution V(a,t) of microscopic 
configurations, we turn to a simplified description, in terms of a finite number of macroscopic 
observabies. The dynamic evolution of these cannot be expected to be closed, as we have lost 
information with respect to the microscopic description. It depends in general on a larger number 
of macroscopic variables, i.e. a hierarchical set of dynamic evolution equations arises. Choosing 
carefully a closure hypothesis at any level of this hierarchy, we are led to improvingly precise 
predictions. 

In the next three subsections the first three levels of this hierarchy are presented, together with 
the corresponding closure assumptions. The evaluation of their accuracy is deferred until Sec.^J 
where we compare them with numerical simulations. 



3.2 The binomial approximation 

The simplest implementation of this idea consists in keeping track of the energy density e(i) and 
of the magnetization per spin m(t) only. 

It is rather natural to include the energy in our set of observabies. Indeed, the system is evolving 
towards equilibrium with respect to the Gibbs measure, in other words, at long times all the 
microscopic configurations with equilibrium energy are equiprobable. Including the magnetization 
is also necessary for a ferromagnetic system, where the low temperature phase is characterized by 
a non zero value of the magnetization. In a more general setting, the minimal set of observabies 
is given by the energy of the system and a complete set of order parameters which allow for the 
exact description of the equilibrium distribution. 

At each time-step, the chosen spin has value a and u unsatisfied edges around it with a certain 
probability p a {u; t). If the flip is accepted, i.e. with probability W(u,(3), the total energy changes 
by an amount of K — 2u (unsatisfied edges become satisfied, and vice versa), and the variation of 
the total magnetization is — 2a. We thus obtain in the thermodynamic limit: 

d K 

— = ^W(u,/3)(if-2u)[p_( U ;i)+P+(w;t)] , (17) 
^ = 2Y,W(u,P)\p-(u;t)-p + (u;t)}. (18) 

u—0 

Although exact, these equations are not of direct use because they involve p a (u;t), a dynamical 
quantity not present in our original description via {e(t),m(i)}. We thus have to express approxi- 
mately p a [u\ t) in terms of e and m in order to close the set of equations. The randomly selected 
variable has spin a with probability (1 +cm(i))/2. In the absence of further informations, we can 
only assume that each of the K incident edges is unsatisfied with the same probability a a (t). This 
yields the approximate binomial expression 

P Au;t) c 1 + g 2 m(t) I* y„(tr(i - a a (t)) K - u . (19) 

ol± can be determined by the following consistency condition. All unsatisfied edges connect anti- 
parallel spins. The energy must then be the same if expressed as the number of + spins around — 
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spins, or the other way around. 

e(t) 



E, . 1 + m(t) Tr . . 



u=0 
K 



Y^up-i^t)^ 1 ™ {t) Ka_{t). (20) 



u=0 



We finally obtain 



1 + <jm(t) /K\ ( 2e{t) \ " / _ 2eft) 



M»?*) = — 2— U(1 + ^)) J " g(i + amft)) J (21) 

This expression has the same form as the equilibrium equation (|llf) . but with e(t) and m(t) being 
dynamical variables, which can differ from their equilibrium values. 
Equations JUJl, l|18() and (|21|l can be condensed into 

de 

- ee F e (e(t)Mt),P) , (22) 

ft 711 

— ee F m (e(t),m(t),(3) . (23) 

A few expected properties of these equations can be checked immediately: 

• F m (e,m — 0,(3) =0 for all e and (3. If the system is strictly unmagnetized, the dynamics 
does not break this symmetry in the thermodynamic limit. In the low temperature phase, 
there is of course an instability with respect to magnetization fluctuations (which are present 
in any finite system), as we shall show in Sec. [S] 

• F e (e p , m = 0, (3) = where e p = 4^(1 + e' 3 ) -1 is the paramagnetic energy density at inverse 
temperature (3. To prove this, one can note that in that case 

p+(K -u) = e- p[K - 2u) p^(u) . (24) 

Using the detailed balance condition on W and the change of variables u h> Jf - u in one 
of the sums defining F e , one thus shows that the paramagnetic state is a fixed point of the 
dynamics. 

• F e (et,mf, (3) — F m (et,mt, (3) = for (3 > (3 C1 where m,f and ej are the magnetization and 
the energy density of the ferromagnetic phase. Indeed, relation (|24|l is still valid, as can 
be seen from (|10|) . and the proof follows the same line as in the paramagnetic case: the 
ferromagnetic state, when present, is also a fixed point of the dynamics. 



3.3 Independent-neighbor approximation 

In the last section, we have described the dynamics of our ferromagnetic model by mapping it to 
two coupled differential equations for the energy and the magnetization. In order to approximately 
close these equations we had to assume a specific binomial form for the quantities p a (u;t) which 
is exactly valid only in thermal equilibrium, but not for intermediate times. It seems therefore 
natural to extend the set of considered observables, and to look for dynamical equations for p a {u; t) 
itself. Note that on regular graphs, these quantities still form a finite set of observables. This would 
not be the case for a graph with unbounded fluctuations of the connectivities. 

To formulate the dynamical equations, we have to take into account different contributions to 
the variation of p,y{u;t) during a time-step: 

• The first type of contribution is due to the flipped spin itself, i.e. a is flipped to — a. All 
incident edges change from satisfied to unsatisfied and vice versa, i.e. we have u <-> K — u. 
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• The second type comes from the neighbors of the flipped spin. For these vertices, a remains 
obviously unchanged, but u is increased or decreased by one depending on whether the 
connecting edge was satisfied or not before the flip. 



Including the corresponding loss and gain terms, we find the following exact equations, 
d 

It 1 



^-Pa(u;t) = -W{u,p) PtT (u;t)+W{K-u,0)p- a (K-u;t) 



+ / ](K - u) W(u,/3) p a (u;t) [— p(u\a, u — y a;t) + p(u - l\a,u — > cr;t)} 

u 

+ u W(u, /3) p-g (u; t) [— p(u\ - a, u — ► a;t) +p(u+ 1| - a,u — ► a;t)] (25) 

u 

where p(u\a,u — > a;t) is the conditional probability that a vertex of spin a belongs exactly to 
u unsatisfied edges, under the condition that this vertex is reached via an edge coming from a 
vertex with spin a and with u unsatisfied bonds. This probability can be computed from the joint 
probability p^ a (u, it; t) of two adjacent vertices which, in the equilibrium context, was already 
introduced at the end of Sec. |21 

p{u\cr,u -> <r;t) = — -. (26) 

As in the previous subsection, these exact equations do not close. The time evolution of the prob- 
abilities p a (u; t) is given in terms of the correlations between neighboring vertices, Pd- a (u-, it; t). We 
can, however, close the equations at least approximately by considering neighbors as independent 



p(u\a, u — y a; t) 



(K - u) p a (u;t) 
(K - u) a 



p(u\-<r,u-KT',t) ~ — , (27) 

i.e. the conditional probability does not depend on the properties of the initial vertex, but only on 
the fact that we reach the new vertex via a satisfied or an unsatisfied edge. In the last equation 
we have used the short-hand notation («) CT = ^ u • p a (u;t). Note that this does not describe an 
average since, for instance, (l)^ equals the fraction of spins with value ct, and normalization holds 
only for (•)+ + (•}_. 

Under this approximation, Eq. Q25J) closes in p a (u;t) and becomes 

^-p a (u;t) = -W{u,(3)p„{u]t) + W{K-u,p)p- a {K-u;t) 
at 

+ ((K- v)W{u^)) a + _ u + 

{K - u) a 

, (uW(u,0)). 



(U)c 



[-up a (u;t) + (u+l) jv(u+l;t)] (28) 



An important observation is that Eq. (|27|l becomes exact for the equilibrium distributions 
(|1(J|) and i|12|) . Therefore the true thermal equilibrium is a fixed point of the closed dynamical 
equations. For intermediate times, however, there will be deviations from Eq. (|28|l . Since the 
description via Pa(u; t) is, however, more detailed than the one of the binomial approximation, we 
expect deviations to be less important. A thorough comparison with Monte-Carlo simulations will 
be given in Sec. 0] 

These equations are ordinary differential equations and can thus be solved numerically by 
standard methods. One has, however, to be careful with the choice of initial conditions. In fact, 
not every normalized p a {u) corresponds to microscopic configurations (<ti, <7n)- A necessary 
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and sufficient consistency condition is given by the fact that each unsatisfied edge connects two 
antiparallel spins, see the discussion before Eq. (|20[1 . Therefore 

Y,up+(u-t) =53«p_(«;t) (29) 

u u 

must hold for arbitrary time t. Restricting the allowed initial conditions to all p a (u;t) fulfilling 
this condition, consistency is preserved by the dynamical evolution l|28[) . This in fact guarantees 
that the only stationary points of the dynamics are the solutions of the equilibrium study of 
Sec. El In the high-temperature phase, the paramagnetic solution attracts the dynamics, while 
the low-temperature phase allows for three stationary points. The two ferromagnetic ones are 
stable, whereas the paramagnetic solution is unstable with respect to any spin-flip asymmetry 
(pa-(u;to) ^ p—cr{u;to)) m the initial condition. 

3.4 Inclusion of neighbor correlations 

In order to further refine the description of the non-equilibrium behavior of the model, we may 
"iterate" the step which led us from the binomial to the independent-neighbor approximation. The 
exact equations for the evolution of the single-site quantity p a (u: t) depend on the joint distribution 
Paa(u,u;t) for neighboring sites. Formulating an equation for the time evolution of the latter 
quantity, which again includes a finite set of functions, we find an equation including three-spin 
correlations. Approximately, these can be expressed in terms of the two-spin distribution, and the 
dynamical equations close. In the formulation of the equations we have to include the effects of 
the flipped spin itself, and of its first and second neighbors. We only give the resulting equations, 
where the time dependence is not stated explicitly to lighten notations: 

■^Pacr(u 1 ,U 2 ) = -W(ui,0) p<r<r(ui,U 2 ) + W(K - «l,/3) P-aa(K - U X ,U 2 + 1) 

+ ^ W(U, j3) [ -p aa {u,Ui) {K -U\- 1) p{u 2 \<7, Ml -> a) 
u 

+p aa (u,u 1 - 1) (K - tii) p(u 2 \cr, ui - 1 — > cr) 
-p_ CT(T (u,wi) (K - ui) p(u 2 \a,ui -> a) 
+P- aa (u,u 1 + 1) [K — u\ — 1) p{u 2 \a,ui +1 — > cr) ] 

+ {ui <-> u 2 ) , 

—p- aa (ui,U2) = -[W(uuP) + W(u 2 ,/3)]p^ aa (u 1 ,u 2 ) 

+W(K- Ul ,P) p aa {K -ux,u 2 - 1) + W(K -u 2 ,/3) p- a -a(ui - l,K-u 2 ) 
+ W ( U >P) [- Po-a{u,Ui) (ui - 1) p(u 2 \ -cr,Ui — > a) 

u 

+p a - a (u,u 1 + 1) u\ p{u 2 \ —a,ui + 1 — ► cr) 

-p-a-a{u,Ui) Ui p(u 2 \ - (T,Ul -> <t) 
+P-a-a{u i Ui - 1) (Ui - 1) p(u 2 \ - CT, Ul - 1 — > u) 

-Paa(u,u 2 ) u 2 p(ui\a, U 2 — ► —a) 

+p aa {u,u 2 - 1) (u 2 - 1) p{ui\a,u 2 -1 -> -a) 

-p- aa {u,u 2 ) (u 2 - 1) p(m\a,u 2 -> -cr) 

+p- aa (u,u 2 + 1) u 2 p(ui\a,u 2 + 1 — > -cr) ] , (30) 

where (ui <-> u 2 ) means that the complete expression on the right-hand side of the equation with u\ 
and u 2 interchanged has to be added. In the following, these equations will be denoted shortly as 
link approximation. As in the case of the single-site quantities, they have to be solved numerically. 
The p a {u) can be recovered from the pair distribution in two distinct ways: 

y,jW ("1,1*2) oc (if - ui) Po-(tti) (31) 

U 2 
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and 

^JV_ CT (ui,u 2 ) cx m p a {ui) . (32) 

These two procedures must yield the same value of p a (u). This consistency condition ensures, in 
analogy to Eq. Ij29(l for the previous approximation scheme, that the only stationary points are 
those given by the static approach. 

In principle, one could go on like this, i.e. equations for higher-order correlations can be written 
down exactly. They will depend on even higher correlations, and an infinite hierarchy of exact 
equations arises. This hierarchy can be cut approximately at any arbitrary level by factorizing 
higher correlations. Since the number of order parameters used to describe the dynamics increases 
with each level, and the observables of lower levels are contained via consistency conditions, the 
description is expected to become more and more precise. On the other hand, the equations become 
more and more complex, as is already clear when one compares them for the three first levels of the 
approximation hierarchy. Instead of continuing further in this way, we compare now the results of 
the first three approximation levels with Monte Carlo simulations. As we shall see, already these 
approximations give an astonishing coincidence with numerical data. 

4 Comparison with numerics 

The approximation schemes presented in the last section do not possess any intrinsic criterion to 
measure their quality. We therefore have performed numerical simulations for large graphs in order 
to check all three schemes. 

To do so, we have first generated large random regular graphs (N = 3 • 10 6 ) of various con- 
nectivities, according to the algorithm of [321 - Then we have performed Monte Carlo simulations 
of the ferromagnetic Ising model defined on this graph. In these simulations, we have used as 
well Glauber as Metropolis dynamics, cf. Sec. 13.11 In the following presentation we concentrate, 
however, solely on the Metropolis case. The results for Glauber dynamics differ quantitatively, 
but the coincidence between numerical and analytical results has the same quality. To suppress 
finite-size fluctuations, which become important for large times close to the critical point, we have 
averaged numerical data for up to 200 independent runs on independently generated graphs. 

In general, we have initialized the system in the following way: all spins are assigned randomly 
and independently a value, with a certain bias in order to impose some initial magnetization on the 
system. Note that this configuration is not an equilibrium configuration for any arbitrary temper- 
ature. Then we have measured the relaxational dynamics from this non-equilibrium configuration 
to thermal equilibrium, and recorded in particular the time evolution of the energy density e(t) 
and of the global magnetization m(t). 

As a result, as represented in Figs. Hand|21 we have found the following behavior: for very 
short and very long time, all three approximation schemes are in extremely good coincidence with 
numerical data. This is even true close to the critical temperature, all three schemes reproduce 
with high precision the critical slowing down, i.e. the diverging longest time scale in the system. 
Exactly at the critical point, all three approximation schemes show the correct algebraic relaxation 
towards equilibrium. So even the simplest approximation, which can be analyzed in large detail, 
see the next section, allows for a precise estimate of the relaxational dynamics close to equilibrium. 

As shown in the insets of Figs. ^ and [21 considerable deviations between numerical data and 
analytical predictions appear only for intermediate times. As to be expected, these deviations 
become smaller for more detailed approximations. In the link-approximation, these deviations are 
hardly visible for all times, even at the critical point. 

In addition, we have observed an increasing precision of all three approximations for growing 
K. We expect therefore, that the binomial approximation becomes exact in the large- if limit |27| . 
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Figure 1: Energy density as a function of time. The three approximations (dashed line = bino- 
mial, dash-dotted = independent neighbor, full = link approximation) are compared to numerical 
simulations (symbols, N = 3 • 10 6 , K = 3, averaged over 200 runs, error bars are smaller than 
the symbol size) for inverse temperatures (3 = 1.0, In 3, 1.2 (top to bottom), with critical value 
(3 C = In 3. In the initial condition all spins were drawn independently, with average magnetization 
0.1. In the main plot, the independent neighbor approximation is not shown because, on the scale 
of the figure, it is very close to the link approximation. The shoulder in the curve for f3 = 1.2 is 
located close to the paramagnetic energy value, the relaxation to the ferromagnetic state appears 
on a longer time scale. The inset enlarges the region of largest deviation between the different 
approximations for (3 C . Obviously, the more involved schemes lead to better approximations, the 
difference between the link approximation and the numerical data is hardly visible. 




t 



Figure 2: Magnetization as a function of time. The parameters and symbols are the same as in 
Fig. H Again, the coincidence of the link approximation with the numerical data is excellent for 
the full time interval. 
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5 Critical behavior 



In this section we investigate in more detail the predictions of our approximations in the critical 
region of temperatures separating the paramagnetic and the ferromagnetic phase. For the sake 
of simplicity, we concentrate on the first two levels of approximations, for which more detailed 
computations can be done explicitly. 

As we have seen before, all the presented dynamical approximations have fixed points corre- 
sponding to paramagnetic and ferromagnetic equilibrium. One has to study now the nature of the 
dynamical flows in the space of projected order parameters. At high temperatures, the only fixed 
point is the paramagnetic one, and it is obvious on physical grounds that it will be stable. In the 
low temperature phase, this fixed point will become unstable and the two ferromagnetic ones will 
attract the dynamics. 

Quite generally, if the projected evolution is described by an n-component vector q(t) of ob- 
servables, the system's evolution is given by n equations <ji = Fi(q) . Fixed points correspond to 
Fi(<lo) = Vt. They are locally stable if and only if all eigenvalues of the n x n matrix A4 are 
negative, where M. is defined by its entries M.y — (dFi)/(dqj), evaluated at the fixed point under 
consideration. In this case, the asymptotic relaxation of a fluctuation towards equilibrium is given 
by the largest eigenvalue X ma x (smallest in absolute value), and the longest relaxation-time scale 
of the model reads 

T = -(Xmax) 1 • (33) 

Approaching the critical point, this eigenvalue tends to zero. The system slows down until the 
relaxation time eventually diverges at p c . Right at the critical point, fluctuations decay only 
algebraically with time. 



dF e 




de 


dm 


dF m 


dF m 


de 


dm 



5.1 Relaxation time and critical slowing down in the binomial approxi- 
mation 

In the simplest case of the binomial approximation, the space of parameters is only two-dimensional, 
the state of the system being characterized by its magnetization and energy. Following the notations 
of Sec. 13.21 we define the matrix M. by 

.11 = | Mee Mem | -= | ■'• ;"" j . c-s-h 

Mme M ran i 

This matrix takes a particularly simple form when computed at the paramagnetic fixed point. 
Indeed, the symmetry of the system under magnetization reversal causes the non-diagonal elements 
to vanish, and one can directly read the eigenvalues of M. in the diagonal entries, 

Mee = ~ +^-f ) K- 1 E (*)w(u,0)e-^(K 2uf , (35) 

' u — ^ ^ 

Mmm = - 2 h *t^ t(*y(u,P)e-^. (36) 



(i + e-n 



It is obvious from these expressions that A4 ee is negative at all temperatures, and that Ai m m 
changes its sign at the critical temperature j3 c = \n(K/(K — 2)). At low temperature, the para- 
magnetic fixed point becomes unstable against fluctuations of the magnetization, as expected on 
physical grounds. 

From the above expression of M. m m one obtains the divergence of the relaxation time of the 
system when approaching the critical temperature from above, 



K - 2 (K 



(37) 
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In the low temperature phase, at the ferromagnetic fixed point (we consider only the one with 
positive magnetization for simplicity), the four elements are different from zero. Both eigenvalues 
of M. are negative, implying the stability of the ferromagnetic phase. As their expressions at 
general temperature is not very illuminating, we concentrate on the (3 — > /3+ limit, for which one 
finds the following scaling of the matrix elements: 



a < , M e 



bVP ~ Pc , M me ~ cy/0 - (3 C , Mmm ~ d((3 - &) , (38) 



where a, 6, c and d can be explicitly computed in terms of W{u, f3 c ). This implies that the relevant 
eigenvalue vanishes as 

\~(d--)(p-p c ). (39) 



After some algebra to simplify the expression, we obtain finally the divergence of the relaxation 
time from the ferromagnetic side of the transition as 



0^02 



(fi - A)" 1 



2{K-2) 



K 



(1 + e -/M" 



(40) 



The two expressions <|37[) and (|40|l show an universal amplitude ratio of 1/2. 



5.2 Relaxation time and critical slowing down in the independent- 
neighbor approximation 

The analysis becomes more involved in the independent-neighbor approximation. As said above, 
we have to compute the 2(K + 1) x 2(K + 1) matrix 



M = 



( dp ai (u^t) 
\dp a2 (u 2 ;t) 



(41) 



starting from equation 1)28(1 ■ evaluated at the equilibrium distribution p a2 { u 2) given in (|1U|) . The 
matrix elements are thus given by 



■M er , -(ui,M2) 



-W(ui,l3) S U1<U2 
{{K-u)W{u,P)) c 



[-(K - m) S UuU2 + (K - tti + 1) S Ul - 



(K - u) a 

+ [-(K - «iK(jh) + (K - til + l)Pa(ui - 1)] 

~(K-u 2 )W(u 2 ,P) {(K - u)W(u, (5))<j (K-u 2 ) 



(K-u) 
(uW(u,i3))- a 



(K-u)l 
—ui <5 M1: m 2 + {u\ + 1) 6 

1il + l,lt2 J 

uW(u,/3))-v 1t 2 



M a -<j{ui,U 2 ) 



+ \-Uip a {ui) + (iti + l)p a {ui + 1)] 

K — U\ ,U2 

, u 2 W{u 2 ,f3) 



(u)l 



[-Uip a (ui) + (Ml + l)p a (ui + 1)] 



(42) 



The relaxation times equal minus the inverse eigenvalues of this matrix. One has, however, to 
be careful since not all eigenvectors correspond to physically allowed deviations of p a {u;i) from 
its equilibrium value. As already discussed in Sec. 13.31 the values of p a {u) are restricted by 
normalization and by the consistency condition l|29|l. The corresponding eigenvectors, one being 
proportional to p a (u) itself, the other one to the deviation dp a (u) / dh according to a deviation of 
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the effective field h away from its self-consistent cavity value, correspond to zero eigenvalues of A4, 
and have to be excluded. 

The most efficient way to achieve this is to explicitly require normalization and consistency by, 
e.g., expressing p±(K;t) through the other values p±(u;t) with u < K. From normalization and 
consistency we immediately find 



Pa{K-t) 



K-l 



K + u K — u 

K P<?W t) + — p-a(u; t) 



(43) 



The matrix A4 becomes thereby reduced to a 2if-dimensional matrix M. with entries («i 2 
0, 11 



dp ai (u!;t) 
dp a2 (u 2 ;t) 



E 



dp^iu^t) dp a (K;t) 



,{ui,u 2 ) - 



t± x dp a (K;t) dp a2 (u 2 ;t) 

K — UG\0 



cr=±l 



2K 



M aua (ux,K) 



(44) 



Besides the excluded unphysical eigenvectors, this matrix has the same eigenvalues and eigenvectors 
as M. 

Unfortunately, we were not able to find a general expression for the smallest eigenvalue of M. 
for arbitrary temperature, and thus of the longest relaxation time of the system. For small values 
of K, the latter can, however, easily be evaluated using a standard computer-algebra system. We 
find, in complete accordance with the binomial approximation, that the longest relaxation time 
scale diverges like r p = Ak((3 c — P)" 1 if the critical point is reached from the paramagnetic side, 
and like ry = ^Ak(P — Pc)^ 1 coming from lower temperature. We thus find the same critical 
exponent and the same universal amplitude ratio 1/2. Only the prefactors are slightly modified, 
with a difference of about 1-2% only between the estimates of the binomial and of the independent- 
neighbor approximations. The values for Metropolis dynamics are recorded in the table below. 



K 


Ak (binomial) 


Ak (independent neighbor) 


3 


tt = 3.2 


^ ~ 3.273043 


4 


§ ~ 0.964286 


ilfgl^ 0.986521 


5 


0.605604 


0.609906 


6 


0.393775 


0.398242 


7 


0.309897 


0.311649 



In Fig. the results of both approximations are compared to numerical data. The results of 
the independent neighbor approximation are slightly higher than the binomial ones. In perfect 
agreement with this observation, also the numerical data are found to systematically deviate to- 
wards slightly higher relaxation times. Note, however, that the extraction of the relaxation time 
in the ferromagnetic phase is slightly subtle: the time interval between the pre-asymptotic and the 
fluctuation-dominated dynamics is rather short even for large systems (N — 3 • 10 6 , averaged over 
20 samples). 



5.3 Algebraic relaxation at criticality 

At the critical point (3 C , the longest time scale diverges, and the observables of the system decay 
only algebraically towards their equilibrium values. Within the binomial approximation, the dy- 
namical exponents for energy and magnetization decay as well as their prefactors can be determined 
analytically. 

The equilibrium at the critical point is given by a non-magnetized state of energy density 
e c = ^k-i) • ^ tb- e initial condition has vanishing magnetization, the energy relaxes exponentially 
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IP-Pel 



Figure 3: Metropolis relaxation time for K = 3 in the paramagnetic (top) and the ferromagnetic 
(bottom) phase. We have plotted the results of the binomial (dashed line) and the independent- 
neighbor (full line) approximations as extracted from the exact eigenvalues of the matrices A4. The 
dotted line gives the asymptotic algebraic divergence in the paramagnetic phase. The symbols are 
extracted from numerical simulations (N = 3 • 10 6 , 20 samples). 



towards e c . However, if the evolution starts from an initial condition of arbitrarily small, but non- 
zero magnetization, the dynamical evolution shows a power-law dependence in time. 

Let us denote the excess-energy density by e = e — e c , and expand the evolution equations 
around their fixed point (e c ,m = 0). Exploiting the spin-flip symmetry m <-> — m and the fact 
that d m F m vanishes at the critical temperature, we obtain for the lowest orders 



de 
~dt 



9 = C§e + C$m i + C$e i + C$em* + ... , (45) 



^ = C{?hm + C { fm 3 + Ci? ) e 2 m + Cifem 3 + C { fm 5 + ..., (46) 
where the coefficients are given by 

1 iV 'l- (, m -M 

(47) 



c (m/e) = J_ d'+ J F m/e (e,TO,/? c ) 

ij i\ j\ <)', dim 



=e c , m=0 



To extract the leading algebraic long-time behavior 

m(t) ~ m t~ Zm , e — e t~ z " , (48) 

we have to compare the dominant terms on both sides of the equations. Seen that the lhs of the first 
equation is of 0(£ -Ze_1 ), the asymptotically larger 0(t~ Zc )-term on the rhs has to be compensated 
by the second contribution with 0(t~ 2Zm ). This results in z e — 2z m . The lhs of the second equation 
is of 0(i~ Zm_1 ), whereas the two dominant terms on its rhs are of 0(t~ z °~ Zm ) = 0(t~ 3Zm ). 
Comparing the order of these terms we consequently find 

z m = ^ , Ze = 1 ■ (49) 

These exponents stand in perfect agreement with numerical simulations. Considering in addition 
the coefficients of the discussed contributions, we are led to 

JF, 1 2 2 F e 
- e °^+2 W ^< (50) 
1 . d 2 F m 1 3 d 3 F m 

-2 mo = eomo »; + 6 m »w' (51) 
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where the derivatives have to be evaluated at the fixed point (e c ,m = 0). After some algebra, we 
obtain the prefactors of magnetization and excess energy: 

1 

(52) 
(53) 

Note that there are two solutions Too = ±^/toq, depending of the sign of the magnetization in the 
initial condition. The quality of the agreement between these predictions and numerical simulations 
(not shown) is comparable to the one of previous subsection. 

6 Two-time correlations 

We briefly sketch in this section an extension of the previously introduced approximations to the 
study of two-time quantities, more precisely of the global auto-correlation function of the spins, 

1 N 

C{t 2 M) = ^$>*(i 2 K(ti) • (54) 

These will be, in the thermodynamic limit, sharply peaked around its average value (with respect 
to the possible histories of the microscopic dynamics). We suppose in the following t 2 > ti to 
simplify the notations. To compute this function, we shall consider t 2 as the evolution time and 
project the dynamics onto a global observable which retains a trace of the microscopic configuration 
at the earlier time t\. More precisely, let us call q ai a-2 (uu, us, su; t\, t 2 ) the fractions of sites whose 
spin equals o\ at time t\ and a 2 at time t 2 , and which have around them 

• uu edges which are unsatisfied at both times, 

• us which are unsatisfied at time t\ and satisfied at time t 2 , 

• su which are satisfied at time t\ and unsatisfied at time t 2 , 

• and ss = K — uu — us — su which are satisfied at both times. 
We introduce the notation 

(•><ri<T2 = E • qa 1 a 2 (uu,us,su;t 1 ,t 2 ) 

uu, us, su 
uu + us + su < K 

which is again not a normalized average for a given value of <j\ and a 2 , but only when the sum is 

taken over all indices, (1)++ + + (1)-+ + (1)-- = 1- The two-time correlation of the spins 

can be obtained from q, and reads in this short-hand notation: 

C(t 2 ,h) = + - - . (56) 

As in the previous cases, one has to impose consistency equations on q, to ensure that the 
number of unsatisfied edges around up or down spins is the same. Enforcing this condition at t\ 
and t 2 yields respectively: 

(uu + us) ++ + (uu + us)^ = (uu + ?is)_ + + (uu + us)__ . (57) 

(uu + su) ++ + (uu + su) h = (uu + su) + _ + (mm + s?i}__ , (58) 
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eo = 



3K (2(K-l) 



2(K-2) 
K -2 



K 



K-l 



K 

E 

L«=0 



u) 



W{u,f3 c )e-^ u 



(55) 




Figure 4: Two-time correlations C(ti + T,ti) as a function of r, for .fT = 3. Solid lines: numerical 
integration of the differential equations, see the text for details. Symbols: Monte Carlo simulations. 
Left: (3 = 1, t\ = 0. Numerical simulations averaged on 200 samples of size N = 10 7 . Right: 
(3 = 1.2, ti = (bottom), 30 (middle), 150 (top). Numerical simulations averaged on 200 samples 
of size N = 3 10 6 , error bars are smaller than the symbol size. 



Boundary conditions also constrain the value of q when t\ = t^. Obviously = q |_ = at 

equal times, and 

q aa (uu, us, su; t, t) = 8 USi o 8 SUi0 p a (uu; t) . (59) 

An evolution equation with respect to ti for q aia2 (uu,us, su;t\,t2) can be closed using factor- 
ization approximations similar to those used in Sec l3.3l here we only state the result: 

d 

— qa ia2 {uu, us, su; t\, t 2 ) = 
dt 2 



— q<j 1 <T 2 { uu i us , su)W{uu + su) + q tTl ,-a 2 ( us i uu , ss)W(ss + us) 

■ [—uu q ai <j 2 ( uu j ms j su ) + i uu + 1)9<ti<t 2 ( uu + 1, ms — 1, su)] 



(uuW(uu + su)) - ai ~<r 2 



{UU) {J1 (j 2 

(usW(uu + stt))_ CT1CT2 



-us q ai a 2 (uu, us, su) + (us + l)q ai<T2 (uu — 1, us + 1, su)] 



us) ( j 1( j 2 

,A\ _ 

su q ai u 2 (uu, us, su) + (su + l)q ax<J2 (uu, us, su+ 1)] 



(suW(uu + su)) ai - a2 



(sU) ai (7 2 

(ssW(uu + su)) ai 



[ ss q ai a 2 (uu, us, su) + (ss + l)q ai(72 (uu, us, su - 1)] , (60) 



where we suppressed the explicit time dependence of q in the right-hand side to lighten notations. 

One can easily exploit these equations to compute the two-time correlation C(t2,t\) for a given 
initial configuration at t = 0. The first task is to determine p a (u;ti) by a numerical integration 
of the equations l|28ll between t — and t = t\. Using then the boundary condition l|59|l and 
the evolution equations H60|) . one obtains q and finally C through (|56|) . We confront in Fig. 21 
the results of such a procedure with Monte Carlo numerical experiments, which are in satisfying 
agreement with each other. 

7 Connection with dynamical replica theory 
7.1 Basic assumptions of dynamical replica theory 

In this section, we are going to show the equivalence of our approximations to the dynamical 
replica theory (DRT) [291 130) under the assumption of replica symmetry. To do so, we will derive 
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the closure equations we used in Sec. starting from the assumptions of the dynamical replica 
approach. 

Similar to our approach, DRT aims at an approximate description of the non-equilibrium dy- 
namics of disordered systems by means of projecting the dynamics onto an array qofn observables, 
and by approximately closing the corresponding equations. The method becomes exact if the fol- 
lowing two properties are fulfilled: 

(i) The observables in q have self-averaging properties along the dynamics, i.e. they follow their 
average trajectories with probability one in the thermodynamic limit. 

(ii) At each time, all microscopic configurations a — (cti, .... , &n) having the same values of all 
observables in q are equi-probable. 

The first assumption is not expected to pose serious problems if, e.g., densities or fractions of 
vertices with certain characteristics are considered. The second assumption, on the other hand, is 
much more crucial: it is obviously true in thermal equilibrium if the energy density is included in q. 
Indeed the probability distribution of microscopic configurations is then given by the Boltzmann- 
Gibbs distribution, i.e. it depends on the actual configuration only via one observable, the energy. 
There is, however, no reason why this should hold far from equilibrium. In fact, as discussed 
before, we have observed that the inclusion of more and more sophisticated observables gives a 
better and better description of the dynamics of the system. 

At variance with the original DRT, where replicas are introduced to average over the disorder 
of the Sherrington-Kirkpatrick model [12], we are going to use the cavity method [7j. The latter 
is more efficient to average over the quenched disorder of diluted systems, in our case represented 
by the random regular graph on which the Ising model is defined. Even if this average is almost 
trivial in the case of a random regular graph, where the disorder enters only via large loops, i.e. 
via self-consistency conditions, this method can be easily extended to fluctuating connectivities or 
random interaction strengths. 



7.2 The binomial approximation 

Let us start the discussion with the binomial approximation, i.e. with the case where the dynamics 
is completely approximated by the behavior of the energy and the magnetization densities. The 
two assumptions of DRT stated above result in the following approximation, 

v(ff t) _ s ( m (t)-jjEi<ri) 8(e(t)-±H(*)) 
[ ' j E* S (m(t) ~jfT,i *<) * (e(t) - ±H(a)) ' (bi j 

Instead of working directly in this generalized micro-canonical framework, we use a generalized 
canonical approach introducing conjugate parameters for the energy, i.e. a formal inverse temper- 
ature (3{e, to), and for the magnetization, i.e. a formal external field j(e, to). We thus replace Eq. 
E3J by 

t) = me J Mejm)) exp | ~P(e, m)H(a) - 7 (e, m) J> } . (62) 

Note that the formal temperature /?(e, to) is different from the physical temperature (3, as long 
as the system is not equilibrated. Both conjugate parameters have to be adjusted such that the 
average values of the energy and magnetization density assume the desired values e(t) and m(t). 
In the thermodynamic limit, the measure becomes sharply concentrated around these values. 

Using these weights for the microscopic configurations, we have to prove that Eq. (|21|l holds, 
i.e. that the two assumptions of DRT lead to the same closure of our approximate dynamical 
equations. This can be easily obtained using the Bethe-Peierls approach sketched in Sec. [3 the 
only modification is due to the additional external field. The analysis follows, however, exactly the 
same steps, and it leads in particular to the desired binomial closure assumption. 
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7.3 The independent-neighbor approximation 



The case of the independent neighbor approximation is only slightly more involved. We have to 
show that the assumption that all configurations with the same p CT (u) are equiprobable leads to the 
desired factorization of the joint distribution of neighbors, and thus to the dynamical equations of 
Sec. 13.31 More precisely, we demonstrate that, under the above-stated assumption, the conditional 
probability of finding a vertex with u 2 unsatisfied edges, following a a — > ±cr edge from a vertex 
with u\ antiparallel neighbors, equals 



p(u 2 | a, ui — > a) 
p(u 2 | a,u\ 



-a) = 



= OC (if - U 2 ) Pa{U2) 

L„ 2 P*A U 1> U 2) 

Pa,-o-(ui,U 2 ) . . 

= (X u 2 p- a {u 2 ) , 

L„ 2 Pa,-a{Ul,U 2 ) 



(63) 



and it is thus independent on u\ , cf . the closure assumptions 1271) . 

Taking any microscopic configuration tr, the distribution p„ (u) is calculated as 



Pcr{u) 



1 N 



(64) 



»=i 



where tij(<j) = £^ • Jij5 crij - crj counts the number of unsatisfied links incident to i. Again, in 
analogy to a microcanonical calculation, we should sum over all configurations having exactly the 
same values of Pa(u) for all a G {±1} and u G {0, K}. As in the the previous section, we can 
circumvent the explicite microcanonical calculation by going to a generalized canonical ensemble 
by introducing formal inverse temperatures (3 a {u) for every pair (a, u). These formal temperatures 
have to be adjusted in order to constrain the p a {u) to the desired values. Out of equilibrium, they 
are not directly related to the physical temperature (3. 

We consequently have to determine the partition function 



N 



ex p\^2^i( u i(^)) \ ■ 

Proceeding in close analogy with Sec. |2 we introduce partial partition functions Z i ij(a. 



(65) 



for 



the subtrees rooted in vertex i, with edges removed. The values of Oi and Uju are fixed, 

where Uju denotes the number of unsatisfied edges including vertex i but not j. Denoting the set 
of all neighbors of i in the subtree by V^u, the partition function can be calculated by an iterative 
procedure: 



icv Mj : \i\=u iH fee/ 



K-l 



^2 Z k \ l ^~<J ll u k \ l ) exp{P-<Ji{u k \i + 1)} 

«fc|i=0 



n 

fcey,| 3 \/ 



K-l 



^2 z k\i(vi,u k \i) exp {Pviiutyi)} 

«ku=0 



(66) 



Introducing generalized cavity fields as hi\j{a, u) = \n{Z i \j{(j 1 u)/Zi\j{—\, K — 1)}, and looking for 
a homogeneous solution h i ^(cr,u) = h(a,u) for all edges we obtain the following 2K — 1 

self-consistency equations: 



l, { n,a) = \n\ X ) ./V.. ,t)" firr.rr)'^ ' " 



(67) 



IS 




h(+,l) h(-,2) 

• • 

+ 



Figure 5: On the left, an example for a subgraph contributing to p_| (2,3) is given. Within the 

cavity calculation, the influence of the exterior part is replaced by the corresponding cavity fields 
h(+,l) and h(-,2), cf. Eq. (JZOJ. 



/(". r) = t—^k— 1 , Z — ■ (68) 



with 

V^" 1 ph(T,«l)+/3-,r(«l+<5<r,-T) 

V K_1 ne' i (+' u i)+/ 3 +("i+ 1 ) 

Note that, by plugging Eq. I|67[) hito l|68|). these equations can easily be reduced to only three 
self-consistent equations for /(+,+), /(+,—) and /(—,—), whereas /(—,+) = 1 by definition. 
Once these quantities are known, we can immediately determine the distribution 



(69) 



As already said, the formal temperatures have to be adjusted in such a way that all Po-(u) take the 
desired values. 

In order to show that the independent-neighbor approximation holds in this ensemble, we have 
to calculate the joint distribution ,0-2(^1,1*2) for neighboring vertices. As an example we are 
analyzing the case a\ = +1, <j<i = — 1, the other combinations work out analogously, see Fig. [5] for 
an illustration: 

P+-( Ul ,u 2 ) oc e M+,«i-i)+M-,«2-i)+/3+(«i)+/M«2) (70) 

\ui - I J \u 2 - I J 

oc m p + {ui) u 2 p^{u 2 )f{+,-Y 1 (71) 

For the conditional probability given in Eq. 1|63|) we thus find 

p(u 2 +, Ui -) = ^ -otu 2 p^(u 2 ) (72) 

which is exactly the independent neighbor approximation applied to close the dynamical equations 
in Sec. 13.31 

This derivation can be easily generalized for the link approximation scheme. We thus conclude 
that the presented approach is equivalent to DRT under the assumption of replica symmetry. 
On the other hand, focusing on finite-connectivity models makes this approach more elegant and 
intuitively understandable, compared to fully connected models for which the analogous equations 
are much more involved 30 . 

8 Conclusions and outlook 

To summarize, we have studied the dynamics of the Ising ferromagnet on a Bethe lattice. This 
simple and, from the point of view of statics, well-analyzed model serves as an ideal testing ground 
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for a series of dynamical approximation schemes first introduced in the context of the analysis 
of stochastic local search optimization algorithms. In particular, we have obtained a detailed 
characterization of the critical behavior of this mean-field model. 

Whereas the presented approximation schemes work very well in this simple case, there remain 
crucial open questions in the dynamic of disordered diluted systems, where the disorder can be 
present either in fluctuating vertex degrees or in randomly chosen interactions strengths. 

A first interesting application of our approach would therefore be the study of disordered 
ferromagnets in their Griffith phase [351 136). This phase is characterized by an anormally slow 
relaxation behavior even in the high-temperature paramagnetic regime, which results from the 
existence of large regions of higher than average coupling strength. 

An even more challenging problem appears if frustration is included into the model, i.e. if 
we turn to systems displaying a low-temperature spin-glass phase. Here the connection of our 
approach with DRT becomes important and opens the way for the inclusion of replica symmetry 
breaking effects. A challenging task in this direction would be to reproduce with such a dynamical 
approach the subtle phenomena of cooling schedule dependence investigated numerically in [37j . 

Finally, we want to mention as a future direction of research the refined analysis of stochastic 
local search algorithms which solve (or approximate) optimization problems like 3-satifiability or 
graph coloring. This includes, e.g., the analysis of the influence of greedy heuristic steps, which 
was out of range before [3"5| . 
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